home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-02 / pnl006.zip / SHCMPLX1.ZIP / SHCMPLX.PAS < prev    next >
Pascal/Delphi Source File  |  1991-03-08  |  10KB  |  357 lines
  1. {$N+,E+}
  2. Unit ShCmplx;
  3. {
  4.                                 ShCmplx
  5.  
  6.        A Complex Arithmetic Unit for Turbo Pascal 5.0 and above.
  7.  
  8.                                    by
  9.  
  10.                    W. G. Madison and Associates, Ltd.
  11.                           8425 Greenbelt Road
  12.                               P.O. Box 898
  13.                           Greenbelt, MD 20770
  14.                              (301)552-7234
  15.                              CIS 73240,342
  16.  
  17.                   Copyright 1991 Madison & Associates
  18.                           All Rights Reserved
  19.  
  20.         This unit and the associated .DOC and TEST*.* files may
  21.         be freely copied and distributed, provided only that no
  22.         fee is charged for the package beyond a nominal copying
  23.         charge, and provided that the package is distributed IN
  24.         UNALTERED FORM. The sole exception to  this  latter re-
  25.         striction is that bona-fide clubs and  user  groups may
  26.         append text material to the documentation file, provid-
  27.         ed that any material appended is clearly  identified as
  28.         to its source, its beginning, and its end.
  29. }
  30.  
  31. interface
  32.  
  33. type
  34.   ComplexElement  = extended;
  35.   ComplexBaseType = record
  36.                       Re,
  37.                       Im  : ComplexElement;
  38.                       end;
  39.   Complex         = ^ComplexBaseType;
  40.  
  41. function CmplxError  : integer;
  42.  
  43. function Cmp2Str(A : Complex; Width, Places : byte) : string;
  44. {Converts a complex value to a string of the form (Re + Im i)}
  45.  
  46. function CmpP2Str(A : Complex; Width, Places : byte) : string;
  47. {Converts a complex in polar form to a string with the angle in radians.}
  48.  
  49. function CmpP2StrD(A : Complex; Width, Places : byte) : string;
  50. {Converts a complex in polar form to a string with the angle in degrees.}
  51.  
  52. procedure CAbs(A  : Complex; var Result : ComplexElement);
  53. function CAbsF(A  : Complex)  : ComplexElement;
  54. {Returns the absolute value of a complex number}
  55.  
  56. procedure CConj(A : Complex; Result : Complex);
  57. function CConjF(A : Complex) : Complex;
  58. {Returns the complex conjugate of a complex number}
  59.  
  60. procedure CAdd(A, B : Complex; Result : Complex);
  61. function CAddF(A, B : Complex) : Complex;
  62. {Returns the sum of two complex numbers A + B}
  63.  
  64. procedure RxC(A : ComplexElement; B : Complex; Result : Complex);
  65. function RxCF(A : ComplexElement; B : Complex) : Complex;
  66. {Returns the complex product of a real and a complex}
  67.  
  68. procedure CSub(A, B : Complex; Result : Complex);
  69. function CSubF(A, B : Complex) : Complex;
  70. {Returns the difference between two complex numbers A - B}
  71.  
  72. procedure CMul(A, B : Complex; Result : Complex);
  73. function CMulF(A, B : Complex) : Complex;
  74. {Returns the product of two complex numbers A * B}
  75.  
  76. procedure CDiv(A, B : Complex; Result : Complex);
  77. function CDivF(A, B : Complex) : Complex;
  78. {Returns the quotient of two complex numbers A / B}
  79.  
  80. procedure C2P(A : Complex; Result : Complex);
  81. function C2PF(A : Complex) : Complex;
  82. {Transforms a complex in cartesian form into polar form}
  83. {The magnitude will be stored in Result^.Re and the angle in Result^.Im}
  84.  
  85. procedure P2C(A : Complex; Result : Complex);
  86. function P2CF(A : Complex) : Complex;
  87. {Transforms a complex in polar form into cartesian form}
  88. {The magnitude is stored in A^.Re and the angle in A^.Im}
  89.  
  90. procedure CpPwrR(A : Complex; R : ComplexElement; Result : Complex);
  91. function CpPwrRF(A : Complex; R : ComplexElement) : Complex;
  92. {Raises a complex (in polar form) to a real power}
  93.  
  94. {===========================================================}
  95.  
  96. implementation
  97.  
  98. const
  99.   MaxStack  = 5;
  100.   StackTop  = MaxStack - 1;
  101.   Pi        = 3.14159265358979;
  102.   PiOver2   = Pi / 2.0;
  103.   TwoPi     = Pi * 2.0;
  104.   F180OverPi= 180.0 / Pi;
  105.   SpConj  : byte = 0;
  106.   SpSum   : byte = 0;
  107.   SpDiff  : byte = 0;
  108.   SpRProd : byte = 0;
  109.   SpProd  : byte = 0;
  110.   SpQuot  : byte = 0;
  111.   SpPolar : byte = 0;
  112.   SpCartsn: byte = 0;
  113.   SpPower : byte = 0;
  114.  
  115. var
  116.   Conj,
  117.   Sum,
  118.   Diff,
  119.   RProd,
  120.   Prod,
  121.   Quot,
  122.   Polar,
  123.   Cartsn,
  124.   Power    : array[0..StackTop] of Complex;
  125.   CmpError : integer;
  126.  
  127. function CmplxError  : integer;
  128.   begin
  129.     CmplxError := CmpError;
  130.     CmpError := 0;
  131.     end;
  132.  
  133. function Real2Str(A : ComplexElement; Width, Places : byte) : string;
  134.   var
  135.     T1  : string;
  136.   begin
  137.     Str(A:Width:Places, T1);
  138.     Real2Str := T1;
  139.     end;
  140.  
  141. function Cmp2Str(A : Complex; Width, Places : byte) : string;
  142.   begin
  143.     if A^.Im >= 0.0 then
  144.       Cmp2Str := '(' + Real2Str(A^.Re, Width, Places) + ' + ' +
  145.                        Real2Str(A^.Im, Width, Places) + 'i)'
  146.     else
  147.       Cmp2Str := '(' + Real2Str(A^.Re, Width, Places) + ' - ' +
  148.                        Real2Str(Abs(A^.Im), Width, Places) + 'i)';
  149.     end;
  150.  
  151. function CmpP2StrD(A : Complex; Width, Places : byte) : string;
  152. {Converts a complex in polar form to a string with the angle in degrees.}
  153.   begin
  154.     CmpP2StrD := '('+Real2Str(A^.Re,Width,Places)+' at '+
  155.       Real2Str((A^.Im*F180OverPi),Width,Places)+#248')';
  156.     end;
  157.  
  158. function CmpP2Str(A : Complex; Width, Places : byte) : string;
  159. {Converts a complex in polar form to a string with the angle in radians.}
  160.   begin
  161.     CmpP2Str := '('+Real2Str(A^.Re,Width,Places)+' at '+
  162.       Real2Str((A^.Im),Width,Places)+' rad)';
  163.     end;
  164.  
  165. procedure IncPtr(var StackPtr : byte);
  166.   begin
  167.     StackPtr := (StackPtr + 1) mod StackTop;
  168.     end;
  169.  
  170. function CAbsF(A  : Complex) : ComplexElement;
  171.   begin
  172.     CmpError := 0;
  173.     CAbsF := sqrt(A^.Re * A^.Re + A^.Im * A^.Im);
  174.     end;
  175.  
  176. procedure CAbs(A  : Complex; var Result : ComplexElement);
  177.   begin
  178.     Result := CAbsF(A);
  179.     end;
  180.  
  181. function CConjF(A : Complex) : Complex;
  182.   begin
  183.     CmpError := 0;
  184.     Conj[SpConj]^.Re := A^.Re;
  185.     Conj[SpConj]^.Im := -A^.Im;
  186.     CConjF := Conj[SpConj];
  187.     IncPtr(SpConj);
  188.     end;
  189.  
  190. procedure CConj(A : Complex; Result : Complex);
  191.   begin
  192.     Result^ := CConjF(A)^;
  193.     end;
  194.  
  195. function CAddF(A, B : Complex) : Complex;
  196.   begin
  197.     CmpError := 0;
  198.     Sum[SpSum]^.Re := A^.Re + B^.Re;
  199.     Sum[SpSum]^.Im := A^.Im + B^.Im;
  200.     CAddF := Sum[SpSum];
  201.     IncPtr(SpSum);
  202.     end;
  203.  
  204. procedure CAdd(A, B : Complex; Result : Complex);
  205.   begin
  206.     Result^ := CAddF(A, B)^;
  207.     end;
  208.  
  209. function RxCF(A : ComplexElement; B : Complex) : Complex;
  210.   begin
  211.     CmpError := 0;
  212.     RProd[SpRProd]^.Re := A * B^.Re;
  213.     RPRod[SpRProd]^.Im := A * B^.Im;
  214.     RxCF := RProd[SpRProd];
  215.     IncPtr(SpRProd);
  216.     end;
  217.  
  218. procedure RxC(A : ComplexElement; B : Complex; Result : Complex);
  219.   begin
  220.     Result^ := RxCF(A, B)^;
  221.     end;
  222.  
  223. function CSubF(A, B : Complex) : Complex;
  224.   begin
  225.     CmpError := 0;
  226.     Diff[SpDiff]^ := CAddF(A, RxCF(-1.0,B))^;
  227.     CSubF := Diff[SpDiff];
  228.     IncPtr(SpDiff);
  229.     end;
  230.  
  231. procedure CSub(A, B : Complex; Result : Complex);
  232.   begin
  233.     Result^ := CSubF(A, B)^;
  234.     end;
  235.  
  236. function CMulF(A, B : Complex) : Complex;
  237.   begin
  238.     CmpError := 0;
  239.     Prod[SpProd]^.Re := A^.Re * B^.Re - A^.Im * B^.Im;
  240.     Prod[SpProd]^.Im := A^.Im * B^.Re + A^.Re * B^.Im;
  241.     CMulF := Prod[SpProd];
  242.     IncPtr(SpProd);
  243.     end;
  244.  
  245. procedure CMul(A, B : Complex; Result : Complex);
  246.   begin
  247.     Result^ := CMulF(A, B)^;
  248.     end;
  249.  
  250. function CDivF(A, B : Complex) : Complex;
  251.   var
  252.     Denom : ComplexElement;
  253.   begin
  254.     CmpError := 0;
  255.     Denom := B^.Re * B^.Re + B^.Im * B^.Im;
  256.     if Denom = 0.0 then begin
  257.       CmpError := -1;
  258.       Quot[SpQuot]^.Re := 0.0;
  259.       Quot[SpQuot]^.Im := 0.0;
  260.       CDivF := Quot[SpQuot];
  261.       IncPtr(SpQuot);
  262.       exit;
  263.       end;
  264.     Quot[SpQuot]^ := RxCF(1.0 / Denom, CMulF(A, CConjF(B)))^;
  265.     CDivF := Quot[SpQuot];
  266.     IncPtr(SpQuot);
  267.     end;
  268.  
  269. procedure CDiv(A, B : Complex; Result : Complex);
  270.   begin
  271.     Result^ := CDivF(A, B)^;
  272.     end;
  273.  
  274. procedure C2P(A : Complex; Result : Complex);
  275. {Transforms a complex in cartesian form into polar form}
  276. {The magnitude will be stored in Result^.Re and the angle in Result^.Im}
  277.   begin
  278.     CmpError := 0;
  279.     Result^.Re := CAbsF(A);
  280.     if abs(A^.Re) > abs(A^.Im) then begin
  281.       {Use the tangent formulation}
  282.       Result^.Im := ArcTan(abs(A^.Im) / abs(A^.Re));
  283.       end
  284.     else begin
  285.       {Use the cotangent formulation}
  286.       Result^.Im := PiOver2 - ArcTan(abs(A^.Re) / abs(A^.Im));
  287.       end;
  288.     if A^.Im > 0 then begin {first or second quadrant}
  289.       if A^.Re < 0.0 then {Second quadrant}
  290.         Result^.Im := Pi - abs(Result^.Im);
  291.       end
  292.     else begin {third or fourth quadrant}
  293.       if A^.Re > 0.0 then {Fourth quadrant}
  294.         Result^.Im := TwoPi - abs(Result^.Im)
  295.       else begin {Third Quadrant}
  296.         Result^.Im := Pi + abs(Result^.Im);
  297.         end;
  298.       end;
  299.     if Result^.Im = TwoPi then Result^.Im := 0.0;
  300.     end;
  301.  
  302. function C2PF(A : Complex) : Complex;
  303.   begin
  304.     C2P(A, Polar[SpPolar]);
  305.     C2PF := Polar[SpPolar];
  306.     IncPtr(SpPolar);
  307.     end;
  308.  
  309. procedure P2C(A : Complex; Result : Complex);
  310. {Transforms a complex in polar form into cartesian form}
  311. {The magnitude is stored in A^.Re and the angle in A^.Im}
  312.   begin
  313.     CmpError := 0;
  314.     Result^.Re := A^.Re * cos(A^.Im);
  315.     Result^.Im := A^.Re * sin(A^.Im);
  316.     end;
  317.  
  318. function P2CF(A : Complex) : Complex;
  319.   begin
  320.     P2C(A, Cartsn[SpCartsn]);
  321.     P2CF := Cartsn[SpCartsn];
  322.     IncPtr(SpCartsn);
  323.     end;
  324.  
  325. function CpPwrRF(A : Complex; R : ComplexElement) : Complex;
  326. {Raises a complex (in polar form) to a real power, returning the result
  327.  also in polar form}
  328.   begin
  329.     CmpError := 0;
  330.     if A^.Re <= 0.0 then begin
  331.       CmpError := -2;
  332.       Power[SpPower]^.Re := 0.0;
  333.       Power[SpPower]^.Im := 0.0;
  334.       exit;
  335.       end;
  336.     Power[SpPower]^.Re := exp(R * ln(A^.Re));
  337.     Power[SpPower]^.Im := R * A^.Im;
  338.     CpPwrRF := Power[SpPower];
  339.     IncPtr(SpPower);
  340.     end;
  341.  
  342. procedure CpPwrR(A : Complex; R : ComplexElement; Result : Complex);
  343.   begin
  344.     Result^ := CpPwrRF(A, R)^;
  345.     end;
  346.  
  347. var
  348.   T1  : byte;
  349.  
  350. begin
  351.   for T1 := 0 to StackTop do begin
  352.     New(Conj[T1]); New(Sum[T1]); New(Diff[T1]);
  353.     New(RProd[T1]); New(Prod[T1]); New(Quot[T1]);
  354.     New(Polar[T1]); New(Cartsn[T1]); New(Power[T1]);
  355.     end;
  356.   end.
  357.